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Effect  of  Nose  Bluntness  on  Twin  Have  Interference  /66 

Xia  Nan  and  Chen  Qiang 
(Department  of  Modern  Mechanics) 

Abstract 

This  paper  discusses  the  mutual  interaction  between  an 
oblique  blast  wave  and  the  shock  wave  layer  on  a  supersonic  blunt 
cone.  A  flow  field  superposition  technique  is  used  to  calculate 
the  shape  of  the  transmitted  blast  wave  and  peak  pressure 
produced  as  it  is  reflected  on  the  cone.  A  detailed  analysis 
shows  that  due  to  the  complexity  of  the  blunt  cone  flow  field , 
especially  the  effect  of  the  high  entropy  layer,  the  shapes  of  the 
transmitted  blast  wave  and  the  peak  pressure  are  quite  different 
from  those  obtained  with  a  slender  one. 

When  a  supersonic  blunt  cone  encounters  a  planar  blast  wave, 
the  planar  blast  wave  meets  the  bow-shaped  separated  shock  wave 
of  the  blunt  cone.  A  transmitted  blast  wave  is  propagating  in 
the  shock  wave  layer  of  the  blunt  cone.  The  transmitted  blast 
wave  is  reflected  from  the  surface  of  the  blunt  cone.  This  is 
the  so-called  "twin  wave  interference". 

The  "twin  wave  interference"  problem  has  an  important 
practical  significance  in  the  aerodynamic  design  of  a  high  speed 
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vehicle.  Since  the  1960's  a  lot  of  work  has  been  done  in 
theoretical  analysis  and  experimental  research.  References  [1] 
and  [2]  give  a  good  overall  summary.  As  for  the  twin  wave 
interference  involving  a  sharp  cone,  because  the  interference 
flow  field  is  similar,  the  problem  is  greatly  simplified.  Good 
results  have  been  obtained  in  engineering  as  well  as  in  numerical 
computation,  in  the  twin  wave  interference  problem  with  a 
supersonic  blunt  object,  however,  because  the  flow  field  of  a 
slender  blunt  cone  is  far  more  complex  than  that  of  a  sharp  cone, 
all  engineering  computation  techniques  are  limited  to  the 
stationary  point  region.  Numerical  calculation is  also 
limited  to  the  blunt  tip  region^r^l.  Moreover,  it  consumes 
quite  a  bit  of  computer  time. 

The  purpose  of  this  paper  is  to  expand  the  flow  field 
superposition  technique  reported  in  the  references  [2,6]  to  the 
interaction  between  an  oblique  blast  wave  and  the  shock  wave 
layer  of  a  blunt  cone.  The  paper  discusses  the  peak  pressure 
variation  on  the  cone  as  affected  by  the  small  blunt  head.  It  is 
also  compared  to  the  situation  of  a  sharp  cone. 

In  order  to  simplify  the  study,  we  only  considered  a  zero 
attack  angle  blunt  cone.  We  also  limited  ourselves  to  the  flow 
within  the  meridian  plane.  It  is  assumed  that  the  gas  is  an 
ideal  gas  with  a  constant  specific  heat,  r*1.4.  All  parameters 
in  the  paper  are  dimensionless.  The  velocity,  pressure  and 
density  are  relative  to  incoming  flow  velocity  ax,  pressure  Poo 
and  density  p^,  respectively.  The  length  is  relative  to  the 
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radius  of  the  spherical  head  r^.  The  entropy  is  relative  to  the 
gas  constant  R. 

I.  Calculation  of  High  Speed  Blunt  Cone  Flow  Field  Parameter 
Before  Interference 


Numerical  solution  is  used  to  represent  the  initial  flow 
field  of  the  blunt  cone  before  interference.  In  the  hybrid  flow 
region  in  the  head  area,  a  linear  method  was  used^7^.  In  the 
rear  supersonic  region,  an  eigenline  method  was  used^8^.  If  x  is 
axial  distance  from  the  body  axis  from  the  center  of  the 
spherical  head  and  r  is  the  distance  from  a  point  in  the  flow 
field  to  the  body  axis,  let  us  define  the  coordinate  €, 


/  6  7 


(1) 


where  Oj&f&l.  The  subscript  s  represents  a  shock  wave  value  and  w 
represents  a  surface  value.  Let  us  transform  the  r,  x  coordinate 
to  {,  x,  and  divide  the  shock  wave  layer  into  N  layers  along  the 
r-direction  and  express  them  as  ^  where  •  •>^n+1"^w“°* 

In  the  x-direction,  let  us  choose  K  straight  lines  where  x  are 
constants  x1<x2<. . .x^.  Then,  let  us  divide  the  entire  shock  wave 
layer  into  meshes  as  shown  in  Figure  1.  We  used  an  intrapolat ion 
technique  to  determine  the  parameters  at  the  nodal  points  (£},Xj) 
of  the  network.  In  calculating  the  supersonic  flow  field  in  the 
rear  of  the  object,  we  only  went  as  far  as  x»15.  Within  the 
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range  of  parameters  calculated,  the  results  agreed  quite  well 
with  other  people's  values^.  Figure  2  shows  the  surface 
pressure  and  Mach  number  distribution  along  the  x-direction  when 
M30b4  and  the  semi-cone  angle  fi>10°.  From  the  figure  we  can  see 
that  the  sharp  cone  has  a  constant  surface  pressure,  while  a 
blunt  cone  has  an  over  exposure  characteristic.  The  lowest 
pressure  is  smaller  than  the  pressure  on  a  sharp  cone.  Then,  it 
gradually  rises  back  to  the  sharp  cone  value.  Because  of  the 
high  temperature  in  the  high  entropy  layer,  the  speed  of  sound  on 
the  surface  of  a  blunt  cone  is  much  larger  than  that  on  a  sharp 
cone.  Hence,  the  Mach  number  of  the  flow  is  much  less  than  the 
sharp  cone  case.  Figure  3  shows  how  the  flow  parameters  vary 
along  the  £  direction  when  x*5.  The  figure  also  indicates  that 
in  the  case  of  a  sharp  cone,  it  is  a  isoentropic  compression 
process  from  the  shock  wave  to  the  surface  of  the  object. 
Therefore,  both  density  and  pressure  are  increased.  As  for  a 
blunt  cone,  because  of  the  high  entropy  layer  on  the  surface,  the 
density  at  the  surface  is  much  lower  than  that  behind  the  shock 
wave. 
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Figure  1  Division  of  the  Flow  Field  Meshes 

1.  planar  blast  wave 

2.  body  shock  wave 

3.  surface 


Figure  2  Pressure  Distribution  on  the  Surface  of  a  Blunt  Cone 
MOC«48*10° 

1.  sharp  cone 

2.  blunt  cone 
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Figure  3  Variation  of  Blunt  Cone  Shock  Have  Layer  Parameters  in 
Three  Directions,  Mgj-4,  8«10°,  x«5 

1.  blunt  cone 

2.  sharp  cone 

II.  Interference  Between  Planar  Blast  Wave  and  Blunt  Cone  Shock 
Wave  Layer 


When  a  planar  blast  wave  traveling  at  a  Mach  number  MB 
encounters  a  supersonic  blunt  cone  at  an  encounter  angle  0B 
(which  is  defined  as  the  angle  between  the  blast  wave  and  a  line 
perpendicular  to  the  body  axis  on  the  meridian  plant  where  v*0) , 
it  is  certain  that  the  blast  wave  collides  with  the  body  shock 
wave  first.  It  creates  two  new  shock  waves  which  are  called 
transmitted  blast  wave  and  transmitted  body  shock  wave, 
respectively.  There  is  a  contact  surface  between  them.  Then, 
the  transmitted  blast  wave  is  reflected  by  the  surface  of  the 
dull  cone  through  the  entire  shock  wave  layer.  Here,  we  will 
primarily  discuss  the  formation  of  this  twin  wave  interference  . 
the  supersonic  region  of  the  cone,  as  shown  in  Figure  4. 


6 


/68 


Figure  4  Encounter  Between  Planar  Blast  Wave  and  Blunt  Cone 

Shock  Wave 

1.  planar  blast  wave 

2.  body  shock  wave 

3.  transmitted  blast  wave 

4.  contact  surface 

5.  transmitted  body  shock  wave 

At  a  certain  instance,  the  planar  blast  wave  meets  the 
curved  body  shock  wave  at  point  Gj.  In  the  body  axis  coordinate 
system  the  velocity  of  Gj  moving  in  the  tangential  direction  of 
the  body  shock  wave  can  be  expressed  as 
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Its  components  in  the  X  and  Y  directions,  respectively,  are 
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where  is  the  angle  of  inclination  of  the  body  shock  wave  at 

Gj,  i.e.  /3.  In  the  body  coordinate  system,  the  parameters  in 
region  C  following  the  planar  blast  wave  can  be  expressed  as 
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The  flow  field  parameters  in  region  1  immediately  after  point  Gj 
have  been  obtained  from  the  numerical  solution.  If  the 
coordinate  system  is  fixed  on  point  Gj,  based  on  the  relation 
that  the  pressure  is  equal  on  either  side  of  the  contact  surface 
and  the  flow  direction  is  consistent,  we  get 


(10) 


..-if  Ml^iayiCoty.-ctgy,  1  .  t,-»[  1 
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where  M^,  A'j  and  M^,  \'c  are  the  Mach  number  and  deflection  angle 
relative  to  a  coordinate  system  fixed  at  for  Region  1  and 
Region  C,  respectively. 


Figure  5  Refraction  and  Reflection  of  Blast  Wave  on  Various 
Interfaces 

1.  blast  wave 

2.  contact  surface 

3.  transmitted  wave 

*1  and  ?c  are  the  angles  of  inclination  of  the  flow  in  region  1 
with  respect  to  the  transmitted  blast  wave  transmitted  body  shock 
wave,  respectively.  Based  on  these  two  equations,  and  ?c  can 
be  determined. 
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With  a  sharp  cone,  the  twin  wave  interference  flow  field  is 
similar  to  the  flow  field  itself.  The  transacted  blast  wave 
Moves  at  a  constant  speed  along  a  line  coving  out  of  the  apex. 

The  shape  of  the  transacted  wave  is  sisiilar  at  different  x 
positions.  In  the  case  of  a  blunt  cone,  the  similarity  does  not 
exist.  At  different  x  positions,  the  shape  of  the  transacted 
blast  wave  is  different.  We  divided  the  entire  shock  wave  layer 
into  N  thin  layers  based  on  the  coordinate  {.  In  a  small  region 
within  to  and  Xj  to  Xj^},  the  transacted  blast  wave  is 

approximated  as  a  straight  line.  The  wavefront  paraaieters  are 
obtained  based  on  the  local  muaerical  solution.  It  moves  m 
parallell  at  a  speed  VGi  (see  Figure  5).  If  a  physical  quantity 
on  the  i**1  layer  is  known,  then  it  is  possible  to  calculate  it  in 
the  i+lth  layer.  If  is  known,  the  angle  between  the 
transmitted  wave  from  the  ith  layer  and  the  x~axis,  ?ix,  is 

(11) 

Based  on  this  equation,  it  is  possible  to  obtain  the  intersect 
between  the  transmitted  blast  wave  and  t}*},  and  the 

coordinates  of  this  point  (ti+i*  Xi<fl) .  Let  ^  be  the  angle 
between  the  tangent  of  ^  and  the  x-axis.  Then,  the  velocity  of 
Gi+1  moving  along  the  tangent  of  *i+1  ,  can  be  written  a8 
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(12) 


.■-i  yi+ft+i  Tt+tn 


(13) 


10 


where 


xjix(i+l)<xj+l 

The  transmitted  blast  wave  is  refracted  after  meeting  the 
interface  at  In  addition,  an  isoentropic  wave  is 

reflected.  There  is  a  contact  surface  between  the  isoentropic 
wave  the  transmitted  wave  layer  in  the  i+l**1  layer.  The  shape  of 
wave  and  the  transmitted  in  the  i+l^  layer  should  meet  the  condition 
that  the  pressure  on  either  side  of  the  contact  surface  is  equal 
and  the  direction  of  velocity  is  consistent.  The  wave  front 
parameters  in  the  ith  and  i+l^*1  regions  are  known.  The  shock 
wave  shape  in  the  ith  layer  is  already  known.  The  parameters 
behind  the  wave  in  region  ia  can  be  determined.  The  shape  of  the 
reflected  wave  and  the  shape  of  the  transmitted  wave  in  the  i+l*-*1 
layer  can  be  found  based  on  two  compatible  conditions;  i.e.  it  is 
possible  to  obtain  To  simplify  the  computation,  a  small 

parameter  expansion  formula  may  be  used  to  derive  an  iterative 
formula.  The  transmitted  waves  of  the  it^>  and  i+lt**  layer  can  be 
obtained  based  on  the  oblique  shock  wave  relation.  The 
isoentropic  wave  employs  a  simple  wave  relationship.  The 
following  functions  are  defined  by  us 


9>), 


(14) 
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The  subscript  a  indicates  parameters  behind  the  wave. 

If  the  coordinate  system  is  fixed  at  the  intersect  Gi+1,  the 
following  equation  can  be  derived  based  on  the  compatibility 
relation  of  the  contact  surface. 


f,  («'<+..  7,9n..)“U-F*(JW*i,y,<p4)  M+i' ^ 

The  parameters  of  two  neighboring  layers  differ  insignificantly. 

Let 


Fl+l-F(+*F, 


(18) 


(19) 


F.(A/,+«.r.9i+ ,)  =  F.(M',,r,9i)+rfF.,  (**1,2,3),  (20) 

By  substituting  these  equations  into  equation  (17) ,  we  get 
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where 
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By  substituting  these  equations  into  equation  (21)  and  changing  it 
from  a  differential  to  a  difference  equation,  we  finally  get 


.  (-h  &y... 

- - - -  OF,  j.  r  _W. 


(22) 


F  i  ”3*.' 
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t*. 


where  K\,  X2  and  1/Pj.  d^\/dn'it  dr2/d*i »  dP2/d*./ 

1/^1  d?i/&p, ,  Fs  are  shown  in  equation  (3-12)  in  reference  (2). 
The  relative  Mach  number  and  inclination  angle  of  the  flow  behind 
the  transmitted  blast  wave  are: 


m',1 


(r-nMy+S  _  + 


2M/’  cot'qp , 

cr- 1)177  +  2 


» 


(23) 
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(24) 


The  angle  between  the  transmitted  blast  wave  and  the  x-axis  is 
expressed  in  equation  (11) . 

As  we  proceed  making  the  calculation  layer  by  layer,  when  we 
reach  the  layer  closest  to  the  cone  surface  which  is  the  N**1  layer 
(see  Figure  6),  the  angle  between  the  transmitted  blast 
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wave  and  the  cone  surface  #a  (which  is  the  incident  angle)  is 
approximately  equal  to  ?N.  if  0  is  the  anqle  of  deflection  of 
the  flow  passing  through  the  transmitted  wave,  it  should  also  be 
the  deflection  angle  of  the  flow  passing  through  the  reflected 
wave. 


(25) 


Figure  6  Reflection  of  Transmitted  Blast  Have  on  the  Surface 
1.  cone  surface 


The  maximum  deflection  angle  of  the  flow  after  passing  through 
the  reflected  shock  wave  is 


(26) 


where  (?,).  is  the  maximum  deflection  angle  of  the  flow. 

The  angle  of  inclination  of  the  corresponding  maximum  reflected 
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shock  wave  is  determined  by  the  following  equation: 


»•"’  v'T  J.- 


r.w. 


Y  -  u-  > 

4"  *"* 


-1  + 


J(y + 1 )(i  +  w;i -r  ItL v#;: )].  (27) 


If  then  it  is  a  routine  reflection.  The  following 

equation  must  be  first  used  to  determine  the  angle  of  inclination 
of  the  reflected  shock  wave. 


-1 


(28) 


The  pressure  behind  the  reflected  shock  wave  is 


if:  > 


(29) 


The  angle  between  the  reflected  shock  wave  and  the  cone  surface 
(angle  of  reflection) 


(30) 


If  it  must  be  a  Mach  type  of  reflection.  The  pressure 

behind  the  wave  is  calculated  by  the  following  equation: 


Pmm  _  tT 

r*  ~yT\*  t+i 


(31) 


There  is  another  case  (Figure  5).  When  the  computation  proceeds 
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to  point  Gi+i#  if  the  Mach  number  is  in  zone  ia  behind  the  wave 
Mja<l  relative  to  the  coordinate  fixed  at  Gi+1#  i.e.  the  flow  ie 
subsonic#  it  is  not  possible  to  produce  a  reflected  wave.  In 
this  situation#  we  can  assume  that  a  Mach  reflection  is  produced 
at  point  G^ .  Make  a  line  perpendicular  to  the  surface  at  G^. 

Its  intersection  at  the  surface  GB  is  used  to  calculate  the 
surface  pressure  behind  the  normal  shock  wave. 

Our  calculation  showed  that  there  is  little  difference 
between  the  results  obtained  based  on  the  iterative  equation  (22) 
and  the  results  obtained  with  a  rigorous  shock  wave  simple  wave 
relationship  by  trial  iterations. 

III.  Calculated  Results  and  Discussion 

Based  on  the  method  discussed  in  this  paper#  we  performed 
computations  on  a  Chinese  made  Model  320  computer  in  two  cases; 
Mao“4 ,  8-10°,  Mb«2  and  M^-5#  8-11.2°,  Mg-2. 
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Figure  7  Oblique  Interference  Pressure  on  a  Sharp  Cone  with  0B 

1.  this  work 

2.  numerical  solution  done  by  Kutler 
W.-5,  M,  — 1.2306,  0-11.2* 


l>t  TIT  38*  «U*  », 


Figure  8  Oblique  Interference  Pressure  on  a  Blunt  Cone  with  0B 
1.  sharp  cone 
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Due  to  the  complexity  of  the  twin  wave  interference 
involving  a  blunt  cone,  we  have  not  yet  seen  any  theoretical  and 
experimental  results  to  date  for  comparison.  For  this  reason, 
the  oblique  interference  pressure  on  a  sharp  cone  calculated  with 
this  method  is  shown  in  Figure  7  together  with  the  numerical 
solution  obtained  by  Kutler^'^J  for  comparison.  The  agreement 
appears  to  be  reasonable.  Figures  8  and  9  show  the  calculated 
peak  pressure  and  geometric  shape  of  the  transmitted  blast  wave 
from  the  interaction  between  an  oblique  blast  wave  and  a  high 
speed  cone,  respectively.  For  comparison,  the  calculated  value 
of  a  sharp  cone  is  also  plotted.  From  Figure  8  we  can  see  that 
the  peak  pressure  on  a  blunt  cone  is  much  lower  than  that  on  a 
sharp  cone  under  similar  conditions.  This  apparently  is  due  to 
the  high  entropy  layer  near  the  surface.  The  high  speed  of  sound 
associated  with  the  high  temperature  makes  the  Mach  number  of  the 
flow  near  the  surface  much  lower  than  that  of  the  transmitted 
blast  wave.  In  the  figure,  X8  is  the  coordinate  of  the  intersect 
between  the  blast  wave  and  the  shock  wave  Gj.  From  the  Figure  we  / 72 

can  also  see  that  when  the  Mach  number  Mw  increases,  this  entropy 
layer  effect  becomes  more  pronounced.  The  peak  pressure  of  a 
sharp  cone  is  independent  of  the  axial  position.  The  flow  field 
of  a  blunt  cone  varies  with  the  axial  position.  Thus,  the  peak 
pressure  also  varies.  The  figure  also  shows  the  angle  0*b,  at 
which  it  changes  from  routine  reflection  to  Mach  reflection.  In 
the  case  of  a  blunt  cone,  it  1b  less  than  that  with  a  sharp  cone. 


18 


Figure  9  Transmitted  Blast  Wave  With  a  Blunt  Cone  in  Oblique 
Interaction  vs.  Body  Axis  Angle  £ 

1.  sharp  cone 

Figure  9  shows  the  geometric  shape  of  the  transmitted  blast 
wave  with  a  blunt  cone  during  oblique  interference.  The  figure 
shows  that  the  slope  of  the  transmitted  blunt  wave  in  the  shock 
wave  layer  of  a  sharp  cone  is  decreasing  uniformly.  In  the  case 
of  a  blunt  cone,  the  shape  of  the  transmitted  blast  wave  also 
varies  quite  a  bit.  This  is  particularly  true  as  the  point  X8 
get 6  closer  to  the  shoulder  of  the  blunt  cone.  As  xs  increases, 
♦x  gradually  decreases  monotonically.  However,  because  the 
entropy  layer  gradient  is  larger  near  the  surface,  ^x  decreases 
faster. 
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When  using  this  method  to  calculate  the  twin  wave 
interference  on  a  blunt  cone,  it  takes  approximately  30  minutes 
to  compute  the  steady  flow  field.  It  only  takes  1-2  minutes  to 
calculate  a  set  of  values  (24  values)  with  a  set  of  given  Xx 
and  MB.  With  respect  to  a  given  set  of  M*  and  6,  the  steady  flow 
field  calculated  does  not  have  to  be  repeated  for  varying  Xg,  mb 
and  6b  in  calculating  the  interference  flow  field. 

IV.  Conclusions 

The  flow  field  superposition  technique  introduced  in 
reference  [2]  was  used  to  solve  the  mutual  interference  between 
an  oblique  blast  wave  and  a  high  speed  blunt  cone  shock  wave 
layer  to  determine  the  geometric  shape  of  the  transmitted  blast 
wave  and  the  peak  pressure  on  the  surface.  The  results  indicate 
that  because  of  the  high  entropy  layer  formed  at  the  high  speed 
blunt  head,  the  entropy  gradient  at  the  surface  is  very  high,  the 
gas  temperature  and  speed  of  sound  are  also  high.  The  Mach 
number  of  the  flow  is  low,  which  makes  the  peak  pressure  much 
lower  as  compared  to  that  with  a  sharp  cone.  In  addition,  due  to 
the  inhomogeneity  of  the  flow  field,  the  shape  of  the  transmitted 
blast  wave  is  also  quite  different  from  that  associated  with  a 
sharp  cone. 
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Application  of  Finite  Element  Method  to  Transonic  Airfoils 
-  A  Preliminary  Steady  Subcritical  Calculation 
Xu  Shoudong  and  An  Changfa 
(Department  of  Modern  Mechanics) 

Abstract 

In  this  paper,  the  finite  element  method  developed  by 
Galerkin  is  used  to  calculate  the  surface  pressure  distribution 
on  a  symmetric  double. curved  airfoil.  Numerical  results 
basically  are  in  agreement  with  experimental  data. 

Based  on  the  transonic  small  perturbation  equations,  the 
finite  element  method  is  used  to  calculate  the  transonic  flow 
field  of  a  symmetric  double  curve  airfoil.  Referring  to 
reference  [1],  we  employed  a  rectangular  mesh  scheme  and  the 
Hermite  interpolation  function.  The  Galerkin  method  was  used  to 
control  the  error.  However,  in  the  finite  element  analysis,  we 
did  not  use  the  numerical  integration  in  determining  the  values 
of  elements  of  the  unity  matrix.  Instead,  we  found  an  iterative 
formula  to  integrate  all  the  elements  in  the  rigidity  matrix. 
Through  the  calculation  of  the  subcritical  pressure  distribution 
of  a  6%  symmetric  airfoil,  it  was  demonstrated  that  the  method  is 
feasible. 

Manuscript  received  on  June  6,  1985 
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I.  Theoretical  Analysis 


The  steady  transonic  flow  around  a  two-dimensional  thin 
airfoil  can  be  described  by  the  transonic  small  perturbation 
equations  and  the  associated  boundary  conditions. 


[l-lg-l&r+l)^]  ^xx+  ^-0 

(1) 

1+  0x)g'*O  (on  the  wing) 

(2) 

^=0,  (at  infinity) 

(3) 

where  g(x)  is  the  dimensionless  wing  surface  equation.  When 
Moo-^l,  it  was  proven  in  reference  [21  that  MooVl+ty+l) 

(local  Mach  number).  Then,  equation  (1)  can  be  written  as 

*xx+*yy-f>  <«> 

As  a  preliminary  approximation,  boundary  condition  (3)  can  be 

considered  valid  at  a  finite  distance. 

If  the  airfoil  is  symmetric  and  the  attack  angle  of  the 
incoming  flow  is  zero,  we  only  have  to  consider  the  upper  half. 

Let  us  choose  a  rectangular  area  which  is  5  times  the  chord  / 75 

length  and  2  times  the  chord  length  height  and  divide  it  into  64 
elements  with  85  nodes  as  shown  in  Figure  1.  For  a  thin  airfoil, 
the  elements  near  the  wing  surface  can  also  be  treated  as 
rectangles.  By  doing  the  following  transformation 
$-(x-xc)/a,  »7«(y-yc)/b,  where  a,  b,  (xc,yc)  are  the  half  length, 
half  height,  and  center  coordinate  of  a  certain  element, 
respectively,  then  all  elements  become  identical  cubic  elements 
as  shown  in  Figure  2. 
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Figure  1 
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Figure  2 

Let's  do  s  Beraite  interpolation  on  the  function  #  in  each 
cubic  eleaent 

(5) 
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where  (i-1,2, . . . ,12)  represent  the  functions  and  derivatives 
on  the  nodes:  •  ^2'  ^2<»  ^2I>'  ^3'  •  ^31* »  ^4 »  » 

♦4,.  M^(^ri})  (i-1,  2,. ..,12)  are  the  shape  functions.  Their 

specific  expressions  can  be  found  in  Conner's  bookie .  For 
example: 


-  07- D(*+l) -(*-!)(*+!)]. 


(6) 


Based  on  the  Galerkin  method,  we  have  the  following  relation 
in  element  E  and  wing  surface  W: 


(7) 


According  to  Green  theorem,  we  get 


\\  &.N,m+i,Nir)d*dy  -j>ti.N,4l-^Mtf..NidBdy 
*  i 


(8) 


where  ?n  is  the  normal  derivative  of  ?  on  the  boundary  of  the 
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element  dE.  The  first  term  on  the  right  side  of  equation  (8) 
happens  to  cancel  out  with  the  first  half  of  the  third  term.  By 
inserting  the  interpolation  function  (5)  and  transformirw  it  into  a 
function  of  ($,*1) ,  we  get 


(9) 


where 


» 


< 


S,,-±An+±B,,,  P, - ±Q 

*  0  r 

f  r  11 

B,,- f  f N,,N,,d;drf 

-l-l  -1-1 

C.I-J  Jy.lV/Mrffrfj;,  X ,-\mN,tdi,  M.+M.+M.+M,), 

-i-i  .  4  . 


(10) 


The  local  finite  element  equation  (9)  is  combined  into  an  overall 
equation 

"is  ,  (U) 

where  k,  1«1,2, . . . ,255.  Then,  we  used  an  iterative  method  to 
solve  equation  (11)  to  obtain  the  nodal  parameters  (including 
derivatives).  Thus,  we  can  find  the  velocity  u«^x,  v-0y  and 
pressure  coefficient  Cp— 20x  at  the  node. 

The  matrices  A^j,  B^j  and  C^j  in  the  element  analysis  were 
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numerically  integrated  on  a  computer  in  reference  [1].  It  will 
bring  about  two  problems:  increased  computer  time  and  new  error 
due  to  approximation.  As  a  matter  of  fact,  all  the  elements  in 
these  matrices  can  be  precisely  integrated.  The  only  problem  is 
there  are  too  many  elements  and  each  element  involves  a  double 
integral.  It  is  just  too  tedious.  However,  an  analysis  shows 
the  following  patterns:  1)  many  elements  in  Gy  are  zero;  2)  all 
elements  of  By  have  a  certain  relation  with  respect  to  the 
elements  of  Ay  and  3)  an  iterative  formula  can  be  found  for  Ay. 
Thus,  the  integration  can  be  greatly  simplified. 

Prom  equation  (6)  we  know  that  the  interpolation  function 
and  its  derivative  are  combinations  of  products  of  (»?-l) ,  (»7+l) , 
(£-1)  and  (£+1)  of  various  orders.  All  elements  of  Ay  are 
double  integrals  of  these  combinations  in  the  cubic  region. 

After  considering  the  separable  variables  of  these  integrals,  it 
is  only  necessary  to  calculate  the  following  integral: 


i>*  (12) 

-I  -I 

By  using  equation  (12),  we  can  calculate  all  the  elements  of  Ay. 
For  example 


{«. 


In  addition,  based  on  the  symaietry  of  ($,*!)  with  respect  to  {, 
V,  we  found  the  following  relations  exist  between  the  elements  of 
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BA j  and  those  of  A^j 


4^8  »  tt  t  -  *•  M  »  (B) 

10  *  II  '‘Vo  0  • 

For  example,  B2,ll“A3,6*  Thus,  when  Ajj  are  calculated,  Bjj  are 
also  determined.  The  computation  of  C^j  is  similar  to  that  of 
A^ j .  it  is  only  simpler.  After  all  the  elements  of  these 
matrices  are  determined,  they  can  be  stored  in  the  computer  for 
future  use  to  avoid  numerical  integration. 

II.  Example  and  Discussion 

In  order  to  check  the  accuracy  and  efficiency  of  the 
technique,  we  used  this  method  to  calculate  two  subcritical 
conditions  of  a  symmetric  airfoil  whose  relative  thickness  is  6% 
at  Ms«0.806  and  0.861.  The  pressure  distributions  thus  obtained 
are  compared  to  the  experimental  data  as  shown  in  Figure  3.  The 
figure  also  shows  the  result  obtained  in  reference  [1].  We  can 
see  that  although  the  meshes  are  divided  more  coarsely  than  in 
reference  [1]  (it  was  divided  into  120  elements  with  150  nodes  in 
reference  [1]),  the  accuracy  of  the  computation  is  similar.  They 
essentially  agree  with  the  experimental  data. 

In  this  paper,  the  computation  is  done  on  a  Model  M-140F 
minicomputer.  Each  state  takes  3  minutes  CPU  time.  In  reference 
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11],  a  mainframe  Onivac  computer  was  used.  Each  state  took  40 
seconds.  Considering  the  difference  in  the  speed  of  these  two 
computers,  the  number  of  operations  of  this  method  is  only  45%  of 
that  used  in  reference  [1].  Thus,  we  proved  that  the  efficiency 
is  improved  with  this  computation  under  the  premise  that  the 
accuracy  is  insured. 

v-  i  ii  ;>  ' 
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1.  this  work 

2.  computation  by  Chan  in  reference  [1] 

3.  experimental  data  by  Knechtel  in  reference  [4] 
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